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1 Introduction 

We consider solutions to the incompressible Navier-Stokes equations 
u t 4- tot* + vu y + wUj + Vp — i/Au, v > 0, 

«,+!!,+«), = 0 , 

on a 2*-periodic square, where u = (u,t>,ti>) is the velocity vector, p the pressure 
and u the kinematic viscosity. Solutions of the Navier-Stokes equations for small 
viscosity are usually iurbuleni ; such flows posses a lot of structure in both space 
and time. The viscosity of the fluid controls the level of turbulence within a flow 
by affecting the energy dissipation. As the viscosity is decreased the size of the 
smallest features, or scales, diminishes. The relation between the viscosity, the 
minimum scale and the total energy dissipation is of fundamental interest for the 
understanding of turbulence. 

The mathematical theory for the Navier-Stokes equations is not complete for 
three dimensional flows: the global regularity is not known and no global bound for 
the velocity gradients is available. However, both these results are known in two 
space dimensions. Whether the results from two-dimensional flows are of physical 
relevance is open to discussion since, in the absence of viscosity, flows in two di- 
mension conserve both energy and enstrophy, while in three dimensions only the 
energy is conserved. Nevertheless, results on two-dimensional turbulence may be of 
significance for large scale oceanographic and atmospheric motions. 

Assuming global regularity, we relate the minimum scale of the flow to iDu^, 
the global bound of the velocity gradients. Our main result, precisely stated in 
theorem (2.1), is that the minimum scale is essentially no smaller than 

*mi» = i ' 1,2 /\D^\ 00 1 ' 7 . 

By comparison, a commonly accepted minimum scale for two dimensional flows, 
^ 2 D i (see Lilly [19], Orszag [20]), is based on the total dissipation rate of the en- 
strophy per unit volume. The enstrophy is defined as the square of the Lj-norm of 
the vorticity. From dimensional arguments it follows that 

a !D = Si'ln'1*, 
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where 

V = 2u J J ||^ || 2 + \\i v \\ 2 dxdy 

is the total rate of enstrophy dissipation per unit volume and ( is the vorticity. 

In three space dimensions the corresponding minimum scale is the Kolmogoroff 
dissipation scale [15] 

A«, = 

where 

e = 2i/ III IK|| 2 + IK II 2 + ||u x || 2 d*dydz, 

is the total rate of energy dissipation per unit volume. 

The estimates for the minimum scale can be used to determine the decay rate 
of the energy spectrum, assuming that a power law does in fact exist. 

From our results in two dimensions we conclude that the energy spectrum, 
E{k ), behaves like k~ s when there is a maximum rate of enstrophy dissipation in 
the flow. The k~ 3 power law is in accordance with the Batchelor-Kraichnan theory 
of enstrophy cascade [3] [16], The high rate of dissipation can not remain for long 
times without the flow disappearing. Indeed, numerical experiments show that the 
solutions rearrange themselves into organised structures which dissipate enstrophy 
at a much smaller rate. Saffman’s work [22], which predicts a power law fe -4 , 
seems to describe the behavior of the system at this later stage of evolution. Our 
theory does not predict the power law but only relates it to the rate of enstrophy 
dissipation; the fc -4 law would correspond to ij of order u 2 ' 2 . 

In three space dimensions there is no a priori bound for ]2>u j^. However, if 
we assume that 

PM* ~ 

then when the energy dissipation rate c is of order one, we obtain the Kolmogoroff 
power law, E(k) = k~ 6 / s , and the Kolmogoroff scale A mtn = A = i/ 3 / 4 . 

Some of the first calculations on two-dimensional turbulence were performed by 
Lilly [19], Fox and Orszag [9], Herring, Orszag, Kraichnan and Fox [12], Fornberg [8], 
and Barker [1], among others. More recent computations on meshes of up to 1024 x 
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1024 points are described in, for example, Brachet, Sulem [7], Brachet, Meneguzzi, 
and Sulem [6], Herring, and McWilliams [11] and Benzi et al [4]. In some cases the 
small viscosity limit of the equations was approximated by the continuous removal 
of the high frequency Fourier coefficients [8]. In other cases the true dissipation 
term is integrated although some extra smoothing of high frequencies is sometimes 
required to suppress the growth of aliasing errors [1]. Another approach is to replace 
the viscosity term by a super-viscosity , that is a higher power of the Laplacian 
operator [7], [11]. This operator allows simulations with a formal viscosity which is 
much smaller. The minimum scale is nevertheless comparable to the computations 
presented in this paper. 

The numerical simulation of three dimensional flows is still limited by the 
power of current computing machines. Currently, the largest three dimensional 
simulations seem to have been performed on 128 s meshes. However, by exploiting 
the symmetries of the Taylor-Green problem, Brachet et al. were able to effectively 
solve with a 256 s resolution [5]. They find the slope of the spectrum to be least 
steep when the rate of energy dissipation reaches a maximum. The numerical results 
seem to agree at this point with the Kolmogoroff scale. For further references on 
three dimensional computations see the review article by Hussaini and Zang [14]. 

We restrict ourselves to two dimensional simulations. Our numerical approach 
has been to attempt to faithfully solve the viscous Navier-Stokes equations. The 
computations were performed using the pseudo-spectral method, Kreiss and Oliger 
[18], and Orszag [21]. There is no extra viscosity added to the numerical simulation 
through smoothing or chopping of the high frequencies, although the fourth-order 
predictor corrector time integrator produces a small amount of it. Numerical simu- 
lations are used to confirm the theoretical estimates and to show that the estimates 
can be achieved for certain initial conditions. Results are shown for the time de- 
velopment of a flow which initially is maximal dissipative. We also show results of 
forced problems. In this case there is no easy a priori bound for the maximum norm 
of the vorticity. We found numerically that the forcing should be proportional to 
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the viscosity in order to obtain order one velocities. More numerical work on this 
subject is still necessary. 

In section 2 we present the analytical results. In section 3 we present numerical 
computations in two space dimensions which substantiate and illustrate various 
features of the proof. Finally, in section 4 we discuss the implications of the 
minimum scale result to the decay rate of the energy spectrum. 


2 Analytical Results 

In this section we will prove some results about the rate of decay of the Fourier 
coefficients for solutions of the incompressible Navier-Stokes equations, 

u t + vu, + vu y + wu z + Vp = i/Au, v > 0, (2.1a) 

u m + v v + w t = 0, (2.1b) 

on the region fi =: {0 < x,y, z < 2ir] and for t > 0. We assume that u = 
(ti(x, <), v(x, t), w(x, *)) is 2x-periodic in x = (e, y, z). 

At t = 0 we give the initial data 

u(x, 0) = Uo(x) , V • uo = 0. 


For simplicity we assume that 


L 


Uo dx = 0, 


which implies 


/„ 


u(x, t)dx = 0 , for t > 0. 


( 2 . 2 ) 


We assume that (2.1) has a bounded solution for all times and want to show 
that the smallest scale is essentially proportional to (iVIjDuI^) 1 / 2 . Here 


|£> U L = sup |£u|oo and |Z?u|oo = sup (|0u/0e|, |0u/0y|, |0u/0z|) . 


In general let 


D p u = 


a p u 


dx Pi dy Pl dz Ps ’ 


I 
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denote any derivative of u of order p, where p = pi + p 2 + ps- 

We are interested in the case when 1 / < 1 and l-Du^ > const > 0. Let us 
further assume that for every natural number p there is a constant C p such that 
the initial conditions satisfy the bounds 


max Hf( 0) =: max ( 

1 0 <j<T 


8’u 

dx> 


0 ; u 


8>u . 


(-. 0)11* + %M)lr) + H^Mr) < c x 


8z> 


Du 

UP 


(2.3) 


Here |f|oo = sup xe0 |f(x)| denotes the maximum norm and 


(f.g)= ft-tdx , ||f|| ! = (f,f), 

J o 

the ^ 2 -scalar product and norm. Then we prove 

Theorem (2.1): We assume (2.3) holds and develop u into a Fourier series 


u(x,t) = £u(k,f)e <<k>x> , k = (fei ,*„*«). 


For every natural number j and any real number a > 0 there are constants K and 
K which depend on j, a and C\ (l = l(j,a)), such that 

y— — rj +1 + 0 


and 


sup |ii(fe,t)| 2 

l>0 


< K l^ u llo +a 

“ \k\ 2 > ' 


The estimate of the theorem can be rewritten in the form 




2j 


We see that the spectrum becomes vanishingly small once |fc| (|Du loo / v ) l/2 with 

|u(fc,<)| decaying faster than any power of (|jDu| 0o /i/) 1 / a |fc| -1 . It is natural to define 
the minimum scale of the flow to be proportional to (» / /|L>u| 00 ) 1/,J . 
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In three space dimensions there is no a priori bound for l-Dul^. One can 
speculate what the right order of magnitude is. For example, if we assume that 

\Dn\oo 

we obtain 

a 

sup|u(fe,<)| 2 < K ^" (^ 3/4 |fe|)~ 2i i 
«>o u 

which corresponds to the Kolmogoroff scale [19] of A m ,-„ = 

In contrast, for two space dimensions an a priori bound for l-Du]^ can be 
obtained. The vorticity, £ = u y — v m , which satisfies 

£t 4- u£ x 4- v£ y = vA£, 


obeys the maximum principle 

\t\oo = sup|£(x,<)| < sup |£(x, 0)j. 

x,t X 

Therefore, assuming that the initial data satisfy 

sup |£(x, 0)1 < 1, (2.4) 

X 

we shall prove that l-Dul^ is essentially bounded independent of i/, in the sense 
that for every J3 > 0 there is a constant K\ = Ki((3) such that 

pu|oo < K x v~ 0 . 

Thus if the initial data have derivatives of order one then the smallest scale is 
essentially of the order i/ 1 / 2 . 

Our proof is also valid for Burger’s equation. In this case one can prove that 
|I?u|oo < const I/ -1 sup |Z>u(as, 0)|oo. 

z 

Thus our result predicts that the minimum scale is of order v~ x . This bound can 
be attained in the presence of shocks. 
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2.1 Estimates for p < 3 

From now on we shall assume that an estimate for exists. Integration 

by parts give us the basic energy estimate 


2 8t 


= -vH 


2 

1 1 


where 


(25) 

Since by assumption (2.3), ||u(-,0)|| 2 < const , it follows 

2v f Hl(t)dt < ||u(-,0)|| 2 < const and ||u(-,f)|| 2 < ||u(-,0)|| 2 < const . (2.6) 
Jo 

Now differentiate (2.1). For any first space derivative Du we obtain 

"II’ + 1 = + l|fu,ll’ + llfu.ll 1 ), 


where 

I = (Du, D(uu x + vu y + tim,)) = II + III , 

II = (Du,itDu a + vDu y + wDu z ), 

III = (Du, Du u„ + Dv Uy + Dw u,). 

Integration by parts and V • u = 0 shows that II = 0. Again by integration by 
parts we obtain 

III < const |Du|ooFf. 


Therefore 


J^llfulp < const | fuU-ff? - •'(llfu.ll* + ||Z>u,||* + ||i>u,|| ! ), 


that is 

~Hl < const iDuUff* - uH\. 

Integrating the last inequality with respect to t gives us, 

H\ (0 < H* (0) + const iDu^ f H\(t)<It — 2v f H\(r)d,T. 

Jo Jo 
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Therefore by (2.3) and (2.6) 

H 2 {t) < const and v J™ H\{t)dt < const (2.7) 

For the second derivatives we obtain 

— + / = -K||D*«.|| J + ||Z> 2 u f || l + ||X> ! u,||’), 

where 

I =(P 2 u, D 2 (uu x + vxij, 4- «m t )) = II + III + IV, 

II =(D 2 u,uD 2 u x + vD 2 u v + wD 2 n z ) = 0, 

III =2(i? 2 u, DuDu x + DvDu v -f DwDu z ) < const |Pu|oo^| > 

IV =(Z? 2 u, D 2 u u a + D 2 v u, 4- D 2 w u z ) < const 

Therefore, 

- const \ Dn \°oHj ~ i 'Hi. 

Integrating the last inequality with respect to t and using (2.3), (2.7) gives us 

2 2 

H\ (0 < cons t — — and v f Hi(t)dt < const -^"j 00 . (2.8) 

v Jo u 2 

For the third derivatives we obtain 

^(Z> 3 u, D 9 u) t + / = -i/(||D s u,|| 2 + ||Z? 3 u, || 2 4- ||Z> 3 u z || 2 ), 
where by Leibniz’s rule 

I = (D 9 u, D 9 {uu x +VU, + wu z )) = 11 + III + IV + V, 

II = (D 9 u, uD 9 u x 4 vD 9 u v 4- wD 9 u z ) = 0, 

III = 3 (D 3 u, Du D 2 u x + Dv D 2 Uy + Dw I? 2 u z ) < const iDu^#!, 

IV = 3 (L> 3 u, D 2 u Du x + D 2 v Du v + D 2 w Du z ) < const iDu^H^Hz, 

V = {D 9 u, D 9 u u* + D 9 v u„ + D 9 w u z ) < const IDuU 
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Therefore 


\ < c»n.( |DuU(H| + - vffj , 

< coni( |Du| w ff| + remit (Zhil^— — ii/if*, 

1/ L 


thus using (2.3), (2.7) and (2.8) 

-s 


Hl(t) < const — ^ - ■ and u J H\(t)dt < const 


(2.9) 


2.2 The estimates for general p 

We now prove theorem (2.1) for arbitrary p. First we obtain energy estimates 
for H p in terms of iD^'ul^, 1 < j < max(l, [(p — l)/3]). These estimates are used 
to obtain bounds for iD^ul^ in terms of IPu^. Finally improved estimates are 
obtained for H p} l-DPu^ and using interpolation inequalities; the theorem 

then follows. We start with 

Lemma (2.1): For every p there is a constant K v such that 
\(D p u>D p (uu m + uu, + mu*)) 

( max(l,[(j»-l)/S]) 

|DnUH, J + J, +1 £ 



Here 


[as] = largest integer < x and l-D^ul*, = max |D*u|. 


Proof: We need to estimate expressions of the form 

(D p VLyD p ~ k u D k u. a + D p ~ k v D k Up + D p ~ h w D k u t ) for k = 0, .. . ,p — 1. 

We integrate by parts to decrease the order of D p u. In doing this the order of 
D p ~ k Uy D p ~ k v and D v ~ k w or the order of D k u», D k n y and D k n z will increase. 
For each new term generated through integration by parts, 

(D q u, D q, u D qs n m ) + (D q u, D q *v D qi u f ) + (D q u, D q, w D q * u z ), 
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we can decrease the order q and increase one of q 2 or g s until one of the following 
conditions is satisfied: 

(1) q - 1 < q 9 < q and q 2 < q 

(2) 9 < 92 < 9 + 1 and 9s < 9 — 1 
Note that in case (1) if 9s = q then 

(D q u, D q, uD q u m + D qa vD q u t + D q *wD q u x ) = 0. 


It follows that 

I = (D p u, D p (un a + vu, -|- wu z ) ) 
can be written as a sum of terms 

A := {D q u, D 2p ' 2q+1 u D q ~ l u, + D 2p ~ 2q+1 v D q ~ 1 u y -f p a r-*«+ 1 tD D q ~ l u M ), 

1 < 2p - 2q + 1 < g, 

B := (D q u, D q+1 v. D 2p ~ 2q ~ 1 u m + D q+1 v D tp ~ Sq ~ l n t + D q+1 w D 2p ~ 2q ~ 1 u z ), 

1 < 2p - 2q - 1 < q - 1, 

C := (D q u, D q u D 2p ~ 2q u. + D q v D 2p ~ 2q u y + D q w D 2p ~ 2q u,), 

0 < 2p — 2q < q — 2. 

Expression -4: If 2p — 2qr + 1 = 1 then the estimate follows, otherwise integration 
by parts is applied to expression A, 

A = (D p ~ q u,D p ~ q+1 (D q uD q - 1 u a ) ) + (D p - q v,D p ~ q+l {D q uD q ' 1 u v ) ) 

+ ( D p ~ q w, D p ' q+1 {D q uD q ~ 1 u s ) ) 

to reduce the order of the factors D 2p ~ 2q+1 it, D 2p ~ 2q+l v and D 2p ~ 2q+1 w. In this 
way we can write A as a sum of terms 

(D q 'u,D p ~ q u + D p ~ q v D qa ~ 1 n y + D*- q w Z? ,,-1 u,), 


where 


p — 9>1, qj < P +1, and 9 i + 92 = p 4- q + 1. 
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Also, 1 < 2p — 2g + 1 < g implies 


p — g < max 



The required estimates can be obtained in the following cases: 

(i) If p — q = 1 then either qi = q 2 = p or one of the qj = p + 1 and the other is 
equal to p — 1. 

(ii) One of the qj is equal to p + 1. 

When neither (i) nor (ii) is satisfied then q j < p + 1 and p — q > 1, and we can 
Teduce D p ~ q further. This shows that A can be estimated in the desired way. 
Expression B: Correspondingly, by reducing the order of D 2p ~ 2q ~ 1 u m , D 2p ~ 2q ~ 1 u v 
and D 2p ~ 2q ~ 1 u z , B can be written as sum of terms 


(D q 'u,D q ’u D p ~ q ~ l n m + D q ’v D p ~ q ~ x u v + D q 'w D p ~ q ~ 1 u z ) := II, 


where 

Qj < P + 1 » 9i + ?2 = P + 9 + 1 and p — q — 1 > 0. 

Also 2p — 2q < q implies 0<p — q — 1 < [|(p — 3)]. II can be estimated in the 
following two cases: 

(i) If p — q - 1 = 0 then qi = p + 1, q 2 = p - 1 or q x = q 2 — p, and 

II < const iDulooHp+iHp-i, or 
II < const \D\x\ooH 2 . 

(ii) If gi = p + 1 or q 2 = P + 1. 

Otherwise qj <p+ 1 and p — q—1 > 0 and hence we can diminish p — q — 1 further. 
Therefore we obtain the desired estimates for B. 

Expression C : Integration by parts allows us to diminish the order of D 2p ~ 2q u 
obtaining terms of the form 

(D qi u, D qa u D p - q ~ 1 u B + D q 'v D p - q ~ 1 u v + D q ’w D p ~ q ~ 1 u z ) 


where 


qj < p + 1 , qi + 92 = P + q + 1 and p - 9 - 1 > 0. 
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Using the same argument as for B we obtain the desired estimate. This proves the 
lemma. 

Differentiating (2.1) p-times gives us 
1 d 

- —(D p u, D v u) + (D p u, D p (u u. + t»u, + wu x ) ) 

= -v(\\D p u m \\> + \\D p u v \\' + \\D p n z \\ 2 ) 

Therefore we obtain from lemma (2.1) 


£ 

;=i 


<amit | IDuUJf, 2 + «,+, Y. I^“l | - 


C mex(l,((j>-l)/S]) \ 

\Da\ x H’ + - Y !*>'“! l«!-i I - "*,’+! 


Using the notation 


the last inequality implies 


= f~ 

Jo 


H l{ 1 ) < const [ 


max(l,[(p-l)/S]) \ 

+ \Dn\„L r +- Y I^IooV/l (210) 


and 


r * I |2?u| |2>uL f 1 

L *+ 1 ^ conat I ~^+r~ + ~^~ Lp + s 


c(l,[(p-l)/S]) 

E 

i = 1 


\ Diu \oo L P -j (2.11) 


To begin with let us obtain estimates for H p and L p+ 1, for p = 4, . . . , 7. In most 
applications this is all what is needed. From our previous results we know that 


La < const 


1 2 H 


and X4 < const 


| 2 ?u| c 


Therefore (2.10) and (2.11) give us 
H\{t) <conat 


Li <conal 


) < 

( | 2 >u|l | 2 ?u| 1 , , 2 \ 

^ L ^ + L_^ X4 + _, I? ul oeZ ; s J 


const 


|Z)u| 


OO 


j 4 ' 


< const 


Ifln.po 


4 

5 " ’ 
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The estimates for H*, L a , H$ and £7 follow in the same way. 

For p — 7 we obtain 

< const ^J^k+lD^,„£r +^15^^+^15*511,^ j- 
Thus we have to estimate |i? 2 u|J 0 . By Sobolev inequalities 

\D p n\l < cH}+ 5 + 


and we obtain for e = H p /H p +s 


\D v n\\ 0 < const H p+ sH p . 


( 2 . 12 ) 


In particular 




< const HeHs < const - — 

Now we use the usual interpolation inequality (see for example [10]) 

l^nlL < e l^ Pu lL + conti e~ (, ~ 1)/(p- - ,) |Z)u|^ ) , 1 < j < p 
which gives us for e = (|£>u|^/|£> p u|^ 0 )( p-, )/( p-1 ) 

|£>'u|^ < const — (2.13) 


For j = 2, p = 3 we obtain 


•13/4 


ID’uli, < const (|i>’a|L) ,/2 (|l>a|L) 1/2 < const (2.14) 


Therefore 


7+1/4 


In | T / 

Hj(t) < const — ^e_ — an d £ g < const 




pi+V* — "■ - „9+l/4 * 

Using (2.14) we can in the same way estimate Hg and if|. These estimates are 
not as sharp as required by theorem (2.1). To obtain the required estimates we 
have to estimate H* for general p and then improve the lower order estimates using 
interpolation inequalities. 
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Using (2.10) and (2.11) recursively and the estimates for H p and £ p+1 for p < 6 
we obtain 
Lemma (2.2): 


ff ; +1 < con., £ 

where the sum is over all {Jt,/,-,p,-,g,} satisfying the constraints 

fc k 


-r* -2A 
OO 


(2.15) 


T.=f+ 1- - 1 ) = 2 * + 


»=1 »=0 

1 < ii < [p</3] 1 Pi+i = Pi - ji - 1 - 9i , 0 < 9i < Pi - ji for i = 1, . . . , k 

0 < 9o < V + 1 , Qk — Pk ~ jh , Pi = P ~ 9o i fc > 0. 

Note that t* = p + 1 /or p < 5 . 

Proof: We first obtain estiamtes for L p+ 2 from (2.11). £T* +1 is bounded by v times 
this estimate for L p + 2 - 

„ J _ 1 •*i**(li[p/S]) 

Z,+j < const ( (Pul^ /^ +J + {\D^\Jv)L p+l + „-2 £ ) 


< const ( 


+ 


B 


+ 


i=l 


). 


where we have labeled the three terms on the right hand side as A, B and C. In the 
recursive reduction of £,42 we must consider all possible terms which may arise; 
at each new stage one must consider the effect of using expression A, B or C. In 
the general case one chooses B g 0 -times followed by term C with j = ji, then B 
gi -times, C with j = ji and so on until finally finishing with term A or B and using 
the estimates for L p with p < 6. In this fashion the general term will be 


\ Diin 'oo 


*1 


\Du\ 


Ik 

OO 


l/lo 




=|0 'hL-I^I’ 


V* 

OO 

00 j.r^+jii +1 


„*» + ! 


|Du| 


We define r* = 2k 4- The constraints on ji, qi and p; follow from the manner 

in which the general term was obtained. 
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Now we use (2.13) to estimate |^' ,< u| 0O in terms of |-D ,, ~ 1 u| 00 and 1-Duj^, 


•®y+i(0 < const ^ |£> p_1 u 

ii—ifc 




r k -2 h 
oo 


V \ | : 2(1 1- i ) 1^11 

< const 2^ |Pi’- 1 u| o ; i 


u Ti 

l+J>(r»-3)/(j>-2) 


i 


)!•••)% 


1/ T * 


Therefore by (2.2) and Sobolev inequalities 


l^-^loo < conat i H l + 1 + ^p-l) 


< const H* +1 


< const max I X» ,_1 u 


=ra*c*-3S?> rDSTL + " T " 5),( '- !) ' 


l/ T * 


It follows that 


|2? p lu loo < const max 


\Du 




r» \ 


One more application of the interpolation inequalities gives for 1 < j < p — 1 


|D J u|„ < const max , , . 

1 lo ° ~ n I i-1+3 ’ 

In order to obtain our final estimate we need to show that that r = min* r* 
tends to infinity as p tends to infinity. Recall that 

fc h 

Tk=P+l~ ]£(;< - !) = 2k + »> 

»=1 i=0 

1 < ji < [p./3] , Pi+i = Pi - ji-l-qi , 0 < $ < pi - ji for « = 1, . . . , k 

0 < 9o < P + 1 , qk = Pk~ jk i Pi = P ~ 9o , * > 0. 

From ji < [p,-/3] < Pi/3 it follows that 

2 

Pi+i > Pi ~ j Pi - 1 - 
(p»+i + 3) > -(?,• + 3) — 9i, 

=* (Pk + 3) > (^) fc - X (p + 3) - 

*=0 
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This last expression can be written as 




Since qh — Ph ~ jh > (2/3 )p* we have 


fc-i 


Tk = 2* + q { > 2k + -p k + ^ ©• 


»=o 

Now consider the two cases 

(i) E,='o( 3 /2)‘?i>(P + 3)/2 

( a ) Eto ( 3 / 2 )'® < 6 > + 3)/2 

In case (i) it follows that 


i=0 


i=0 

*-i 


A-l 

^(3/2) fc - 1 g< > £(3/2)* g,- > (p + 3)/2 

i=0 


=* ^ 9i >(2/3)*- 1 (p + 3)/2 
*=o 

\h-li 


=>r h >2k + (2/3)" -1 (p + 3)/2 
Minimizing this last expression with respect to fc (> 0) gives 

Tk > const log(p) for some const > 0. 


In case (ii) we have 


* * 1 + I °«>/ j( 2(STs) ) 

=> 71 > 2 + 21og a/ ,(j^i/jj) + (2/3)pi 

Since 0 < p* < p + 1 it follows that as a function of p*, the above expression has 
the bound 


t* > const log(p) for some const > 0. 

Hence r = min* t* > log(p) as p tends to infinity, for some constant 6 > 0. Thus 
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Theorem (2.2): For every j > 1 and any a > 0, we can choose p sufficiently large 
so that 


|D>u < const 


|Du 


r ;' + l+Q 


(2.16) 


1/ ) — 1 + <* 

where the constant depends on p and C P+ i introduced in the estimates for the initial 
conditions (2.3) . 

By using (2.16) in (2.15) we obtain 

Theorem (2.3): For every j > 1 and any a > 0, we can choose p sufficiently large 
so that 


|Du 


ri+a 


(2.17) 


Hj(t) < const — T-r — 

3 W — if J+<* 1 

where the constant depends on p and c p+ i . 

Using the simple estimates for Hj in terms of maximum norms gives 
Theorem (2.4): For every j > 1 and any a > 0, we can choose p sufficiently large 
so that 

y— — ri + l+a 

Hf{t) < const . (2.18) 

where the constant depends on p and C P + 1 . 

Theorem (2.1) now follows from Parseval’s relation. 

It may seem curious that the initial conditions satisfy 

27/(0) < const -y 22 -, 

while we are able to prove that (2.18) holds. However, (2.18) can be derived from 
(2.17) as follows (for convenience we drop the a’s) : 

phL 


H; < const 
1 - Vi 


=$■ |D*u|^ < const H* +i < const 


| Du 


rJ>+* 


j >+2 


and thus using the interpolation inequality (2.13) 

|D'u|^ < const 

< const ( 1— — W— )( , ~ 1 )/( p ~ 1 )(|Du|| 0 )( > ~ , )/( y ~ 1 ) 


< const 


| Du 


i/P+ 2 

•J+l+a 

oo 

i — 1 + a 


- 18 - 


(2.18) now follows. 

2.3 Estimates for two space dimensions 

In this section we obtain a sharp bound for l-Du^ for two dimensional flows. 
In this case the incompressible Navier-Stokes equations can be written in vorticity 
form 


it + + vil = vAi , v > 0, (2.19a) 

u m + = 0 , — v m = i, (2.19b) 

where ( is the vorticity and u , v are the velocity components. 

Lemma (2.2): The solutions of (£. 19 ) satisfy the maximum principle 

l<(-.0loo <I*M)U. 

Proof: This well known result follows from the fact that at a local maximum 
(minimum) of i, i* = i v = 0 and & < 0 (> 0). 

We would first like to show that l-Dul,* is bounded for all time. Note that our 
energy estimates of the previous section are still valid if we integrate to t — T > 0 
instead of integrating to t = oo. For this section only, let us redefine the quantities 
which depend on this bound on the time. For example we define 

f T 

iDul^ = suplDuloo and L p = H p (t)dt. 
t<T Jo 

We know from basic results that |Du|eo exists and is bounded for some finite time 
interval [0,T]. We will now derive estimates for l-Du^ which are independent of 
T. It follows from the results of the previous section that we can obtain estimates 
for all derivatives which are also independent of T. Then from well known results 
we can conclude that jDuloo exists and satisfies these same bounds for all times. 
Lemma (2.3): For any ati > 0 there exists a constant C(ati) such that 

1KI. = ■npIBnU < C(«i)l((->0)|! o + *V— • 

t<T 


( 2 . 20 ) 
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Proof: For any /3 with 0 < j9 < 1 we define the Holder semi-norm by 

= 811D l/( x 0-/( x 2)l 

x,,xfeO |*1 -X 2 |^ 

Using the notation 

iZ^uloo = max{|Z? /3 it| 0C) , \D 0 v |oo} 

the usual Holder estimates for the solutions of Laplace’s equation (see for example 
[10] ), tell us that for any /3 > 0 there is a constant C(/3) such that 



|C‘-'uU < C{0)\(\ x 


(2.21) 


Also, the convexity of Holder norms (see [13]) 

|£* +t 7U < const |P* 1+ai /l«|^ a+aa /|«‘ , 

fc + a = t(&l + Qtl) + (1 — t)(k 2 + aj) > 1, 0 < t < 1, CK, cri, CK 2 > 0, 
and Young’s inequality give us for any e > 0 


|0u|oo < clP^loo + cons* e p \D l ^ul*. 


Using the Sobolev inequality (2.12) for |L >2 u| 0o and (2.21) we obtain 

nri 7/2 

l- Du loo ^ conti e + const e - ^|P l_ ^u] 00 , 

un 7 / 2 

< const | e + e _/} C(/3)]^i 00 ) . 


Choosing 


gives 


Mfi _ 


= com 


z 7 / 2 


i^iL /2 ' 



i / 7 / 2 \ 

l^[L /2 / 


“TTF 
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Thus 

\D^L<am.t (dm") 1 ' 11 ''’ 1 " 

and the lemma follows since 1^1^ < |£(’,0)|oo. 

In two space dimensions estimates on the vorticity appear more naturally. In 
[17] we proved the results of theorem (2.1) in the two-dimensional case using the 
vorticity formulation of the equations. In that paper the quantities 


m = 




+ 


.n. 


0 **" "dir" 

take the place of the H*. The estimate corresponding to (2.18) is 

-rJF-f 2 + a 


( 2 . 22 ) 


Jp(t) < const Hp +1 (t) < const 


|Z?u| 




(2.23) 


We refer to the J r in the section on numerical results. 
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3 Numerical Results 

We first describe the procedure we use to numerically solve the two-dimensional 
Navier-Stokes equations. In brief, we discretize in space using the Fourier (pseudo- 
spectral) method and solve in time using a fourth order predictor-corrector method. 
The equations are solved in Fourier space and the diffusion term uAu is treated in 
a fully implicit manner. We now proceed to present more details. 

We solve the two-dimensional incompressible Navier-Stokes equations in the 
vorticity stream function formulation: 

it + («$)■ + (»£)* = v ^i + / ( 31a ) 

A Y> = -i, («,*>) = (3.1b) 

The computational domain is taken to be a 2ir periodic square. The solution is 
represented as a truncated Fourier series with denoting the discrete approximation 
to £ and <2> denoting the discrete Fourier transform of u>: 

»■(*.»,<)= Y, T, *(* 

“ 5 - ATi + 1 “ j + 1 

Similarly the Fourier transform of if> and / are denoted by and / respectively. 
The equation for the Fourier coefficient u>(fei,fe 2 ,t) is 

u >t + *A 1 («u>) + ikiiyu) — — i/(As* -I- kf)u> 4- f (3.2a) 

+ k\)j> = u> (3.2b) 

The convolutions rtu; and vu (i.e. the Fourier transforms of the products uu> and vw) 
are computed from «, v and u> by transforming to real space, forming the products 
and then transforming back to Fourier space, (pseudo-spectral method). It is not 
hard to see that the computation of u>j can be done with five two-dimensional fast 
Fourier transforms (FFT’s). In fact, only four FFT’s are needed since one can write 
(reference Basdevant [2]) 


V’» U, K (V’j) )»» [(^» (^’»V’jr)lf»]* 


- 22 - 


However, in the calculations presented here the less efficient method was used. 

The equations (3.2) can be written in the form of a large system of ordinary 
differential equations: 

F(y,f) (3.3) 

where y is the vector with components b>(&i,fc 2 ). 

Time stepping is performed using a predictor-corrector applied directly to equa- 
tion (3.3). Let y w denote the approximation to y(nAt) and F n = F(y n ,nA<). We 
use the fourth order Adams predictor-corrector scheme given by 

y p = y» + y^(23F n - 16F n _! + 5F n _ 2 ) (3.4a) 

yn +1 = y» + ^( 9F r + 19F » - 5F »-1 + f »- 2 )- (3.4b) 

Here y f is the result of the Adams-Bashforth predictor, F p = F(y F ,(n 4- l)At) 
and y n +i is the corrected value obtained from approximating the implicit Adams- 
Moulton scheme. A single time step thus requires two evaluations of the right hand 
side F. The classical fourth order Runge-Kutta method is used to obtain starting 
values for (3.4). These are required initially and whenever the time step is changed. 

For stability reasons one may want to integrate the diffusion term, i/Aur, in an 
implicit manner. In the Fourier representation this term is very simple and thus can 
be easily treated in a fully implicit and accurate manner. We write the equations 
(3.3) in the following way 

^ = G(y , t) - Ay A = diag(..., v{k\ + k \), ...). 

where the right hand side F has been split with A the diagonal matrix corresponding 
to the diffusion term. This last equation can be written in the form 

^(e A ‘y) = e A *G. 

Now apply the time stepping procedure (3.4) to this equation viewed in terms of the 
new dependent variable e A, y. After division by e A< the predictor-corrector scheme 
which results is 

y„ = e~ AA< y» + ^(23e- AA ‘G n - 16 e - 2AA *G n _ 1 + 5 e - SAA< G n _ 2 ) (3.5a) 
y»+i = e- AA ‘y« + ^(9G„ + 19 e - AAt G n - 5 e - JAA< G n _i + e - 3AAt G n _ a ).(3.5b) 
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The Runge-Kutta scheme is transformed in a similar fashion. The terms e - ^ fc > +fc ^ A< 
are stored and need only be recalculated when At changes. These resulting schemes 
are exact in the absence of the convection terms (G = 0). 

The variable time step is chosen by stability and accuracy considerations with 
At chosen to satisfy the condition 

CFL m in < (|«U + Moo)^ < CFL max (3.6) 

h 

where h = 2/N, (N = max(iV 1 , 7NT 3 )). The stability region of the explicit predictor- 
corrector method (3.4) is shown in figure 1 . 



Figure 1 Stability region for the predictor-corrector scheme. 

When (3.4) is applied to the model problem y' = Ay, the time step is restricted 
by (approximately) |A|Af < 1.2 if A is purely imaginary and by — AAt < 1.9 if A is 
real. One expects the implicit predictor-corrector scheme (3.5) to have better sta- 
bility properties than the explicit one (3.4). CFL m * m and CFL max are the minimum 
and maximum allowable values for the Courant-Friedrichs-Lewy number: 


CFL :=(Moo + Mco)^- 

h 
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CFL max would be taken less than the stability limit for the model problem. (The 
choice of h in our definition of CFL instead of the true h = 2i r/N means that we can 
compare CFL directly to the normal stability limit for the model problem.) When 
the condition (3.6) is violated the new time step is chosen so that 

At 

(Moo + Moo)-f~ = CFL opt . 

h 

3.1 Verification of the Numerical Approximation 

In this section we present results which illustrate the accuracy of the numerical 
approximation that we use. In test 1 we show that the time stepping procedure 
is accurate to fourth order in At. In test 2 we consider the convergence of the 
numerical solution as the number of modes is increased. 

Test 1: Accuracy of the time stepping procedure 

It is easy enough to choose the forcing / in the Navier-Stokes equations (3.1) 
so that the true solution is known to be some given function. Numerous tests of this 
kind were performed. In all cases the numerical solutions converged to the exact 
solutions at a rate very close to fourth order in the time step At. 

As a more realistic study of the convergence of the time stepping routine we 
consider a sequence of calculations with fixed random initial data and decreasing 
time steps. The initial conditions are identical with those used in section (3.2) for 
the decay of random initial data. Keeping the same initial conditions, and with 
JV - ! = Nz = 128, u = 10~ 4 , the equations were solved with three different (fixed) 
time steps: At = .05, At = .025 and At = .0125. The computed maximum value 
for the CFL number in each run was 1.2, .6 and .3, respectively. (Recall that the 
stability limit for the explicit version of the predictor-corrector scheme is about 1.2 
on the imaginary axis. We were able to obtain good results for values of the CFL 
number as large as 1.5, which substantiates the belief that the implicit predictor- 
corrector scheme has better stability properties.) We use the results from the three 
runs to estimate the rate of convergence as a function of At as well as to estimate 
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the actual error. We measure both the discrete maximum error and the l 2 error 
defined by 

a>|oo = max |u> 

(»<.*>) 

By assuming that the computed solution is converging to the true solution as 0(At p ) 
we can determine approximate values for p and the errors 


(*<.»)! and ^ K*. -,%•)! 




Eoo(t, At) := |u; coinputed (-,<;A<)-u> trU e(-,<; A<)|oo = 0(At p ), 
E 2 (t, At) := |a; com p uled (-, <; At) - a> true (-, <; A<)| 2 = 0( At p ). 


These values are given in table I for the maximum norm errors and in table II for 
the l 2 errors. 


t 

Eoo (At 

= .05) 

EooiAt : 

= .025) 

£00 (Ai 

- 

.0125) 

V 

10. 

0.68 

X 

10- 2 

0.47 

X 

10~ 8 

0.32 

X 

10- 4 

3.9 

20. 

0.27 

X 

10~ 2 

0.18 

X 

10" 8 

0.13 

X 

10- 4 

3.9 

30. 

0.97 

X 

io- s 

0.69 

X 

10- 4 

0.49 

X 

10" 8 

3.8 

40. 

0.69 

X 

10" 8 

0.47 

X 

10- 4 

0.32 

X 

10~ 8 

3.9 

50. 

0.48 

X 

10- 8 

0.31 

X 

10- 4 

0.20 

X 

10" 8 

4.0 


Table II - Estimated maximum errors and convergence rate: 0(A< P ) 


t 


= .05) 

E2 ( - 

= .025) 

E 2 (At 

= 

.0125) 


10. 

0.40 

X 

10- 8 

0.29 

X 

10- 4 

0.21 

X 

10- 8 

3.8 

20. 

0.17 

X 

10- 8 

0.12 

X 

10- 4 

0.86 

X 

10" 6 

3.8 

30. 

0.73 

X 

10- 4 

0.50 

X 

10- 5 

0.34 

X 

10" 8 

3.9 

40. 

0.56 

X 

10- 4 

0.37 

X 

10" 8 

0.24 

X 

10" 8 

3.9 

50. 

0.48 

X 

10- 4 

0.31 

X 

10" 8 

0.20 

X 

10-® 

4.0 


Table II - Estimated l 2 errors and convergence rate: 0(At p ) 


Test 2: Convergence for random initial data 

We consider the computation which is described in section (3.2) under the 
heading of Run 1: Decay of Random Initial Data. This computation was run with 
N = N\ = N 2 — 128 (u;i2s) and also with N — 256 (u> 2 5 e)- The initial conditions for 
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the two runs were the same to single precision (about 6 — 7 decimal digits), although 
the actual computations were done in double precision. The variable time step was 
determined by the parameters (CFL m j n , CFL op ti CFL ma x) = (.8,1., 1.2) In table II 
we indicate the maximum difference and the Z 3 difference between the two runs at 
various times. Due to the variable time step the solutions were not compared at 
exactly the same times. The difference between the times is given in the table as 
<258 — < 128 * Note that the maximum difference between the two solutions occurs at 
smaller times. Later on when the solution becomes smoother the errors are smaller. 
Further details of this run can be found in the next section. 


t 

I w 258|<» 

1^258 ~ W 1M oo 

|“>258 ~ w 128 2 

<258 ~ <128 

■ 

0 . 

1.00 

.10 x 10- 8 

.38 X 10“ e 

0.0 

Bfl 

.75 

.55 x 10" 1 

.45 x 10- 2 

+ .42 X lO" 8 

mm 

.67 

.20 x lO' 1 

.29 x lO" 2 

-.21 x 10" 3 

30. 

.60 

.15 x lO" 1 

.21 x lO" 2 

+.23 x 10" 3 

40. 

.59 

.10 x 10~ l 

.73 x 10" 3 

+.18 x 10" 4 

50. 

.58 

.73 x 10“ 2 

.58 x 10~ 3 

+.23 x 10" 4 

60. 

.57 

.43 x 10' 2 

.88 x 10" 3 

-.19 x 10" 3 

70. 

.56 

.51 x 10" 2 

.88 x 10" 3 

+.25 x 10" 3 

80. 

.56 

.26 x 10“ 2 

.38 x 10" 3 

+.37 x lO" 4 

90. 

.56 

.27 x 10" 2 

.36 x 10' 3 

+.12 x 10~ 4 

100. 

.55 

i 

.31 x 10~ 2 

.37 x 10~ 3 

-.62 x 10~ 4 


Table II - Convergence of random initial data v = 10 4 
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3.2 Computational Results 

In this section we present the results of four different runs: 

(1) Run I: Decay of random initial data, v — 10 -4 , N = 256 and N = 128. 

(2) Run II: Decay of random initial data, v — 10 -5 , N = 512. 

(3) Run III: Decay of smooth random initial data, v = 2 x 10 -B , N = 256. 

(4) Run IV: Random forcing, i/=.5x 10 -s , v = .5 x 10 -4 . 

Run I: Decay of random initial data, v — 10~ 4 , N— 256 and N=128. 

For the first run we consider the time evolution of the Navier-Stokes equations 
for random initial data. The initial conditions for the vorticity were chosen so that 

with a random phase. (Actually the initial spectrum was set to zero for all wave 
numbers above some large value of k.) The constant C was determined by nor- 
malizing the maximum value of the vorticity to be 1 at t = 0, |u>(-, •, O)!,*, = 1. 
The value of the viscosity was taken as u = 10 -4 and the number of modes was 
N\ = N 2 = 256. We show 

1) contour plots of the vorticity (figure 2 ). Dashed lines indicate negative contour 
values. 

2) surface plots of u>(Aci, Asj) in the cosine-sine representation. The discrete Fourier 
series for u is actually represented in the computer code as a real series in 
cosines and sines. The surface plot shows the magnitude of the coefficients of 
this series. The coefficients are ordered in the following manner: 


‘Cl Cl 

Cl«l 

C1C2 

Cl *2 

C1C3 ..." 

*lCi 

*1*1 

*1C 2 

*1*2 

*lC S . . . 

CjCi 

Cj*i 

c 2 c 2 

C 2 *2 

CjC 3 . . . 

SjCi 

*2*1 

* 2 c 2 

*2*2 

S2C3 • • • 

C3C1 

Cs*l 

CsC 2 

Cs*2 

C3C3 . . . 

: 

• 

; 

; 



where c*cj is the coefficient of cos(kse) cos(ly), Ck»i the coefficient of cos(kz) sin (ly). 
and so on. The lowest frequency modes are located at the top of the surface 
plot (figure 3 ). Only the first 128 modes are shown in the surface plots. 
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3) plots of the energy, enstrophy, as a function of time, and the decay 

of the vorticity spectrum as a function of k. In figure 4 a we plot the square 
root of the total energy 

(i/2(M ! + imh 1 /*. 

the square root of the enstrophy 


|u>| 


) 


and 

V */V lW = v V»(|| w .||» + | K ||*)«/» 


as functions of time. In figure 4 b we plot the normalized versions of t/ 1 / 2 /i(t), 
t/ 2 / 2 7 2 (<) and u^ 2 Js(t). Recall that 




In each case the functions plotted are scaled so their maximum value is 1. This 
maximum value is indicated on the plot as the value of Scale. In figure 4 c some 
selected Fourier coefficients are plotted as functions of time. Finally in figure 
4 d we show log-log plots of u>(k) versus k. The quantity v(k),k = 1,2,... 
is defined to be the average value of |ui(Z 1 ,/ 2 )| over all wave vectors (fi,I 2 ) for 
which k is the closest integer to l = |(Ii,Is)|: 

*(*)=( E E J ) 


We plot log 10 (<2>(fe)) versus log 10 (A) for different times. Lines with slopes — 1 
and —2 are also marked. Note that if u>(k) ~ k~° then E(k) ~ k~ 2a ~ l . 


For comparison, in figures 5 - 6 , we show the results of the same run when only 
half as many modes were used, N = 128. Essentially the only noticeable difference 
is in the plot of the spectral decay. A quantitative comparison of the 256 and 128 
runs was given in section (3.1). 


Run II: Decay of random initial data, v = 10 -B , N = 512. 
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The value of v is taken as 10 -6 . The initial conditions are the same as Run I. 
The number of modes was Ni = N 2 = 512. The results are shown in figures 7 - 
8 . Note that for technical reasons the contour plots were made by projecting the 
solution to a 256 x 256 grid. It is interesting to compare this run (i/ = 10 -s ) to the 
previous run (i/ = 10~ 4 ). 

Run III: Decay of smooth random initial data, u = 2 x 10 -5 . 

In this run we begin with initial data which is much smoother than in the 
previous runs. The initial vorticity spectrum is chosen so that 

|&(Ai,* 2 )| = Cite" ’ (fc/fco)S , *0 = 3.5 

with random phase. The constant C is chosen so that |u>(-, •, O)),*, = 1. The viscosity 
was 2 x 10 5 and the number of modes was N = 256. These initial conditions are 
similar to those used by Brachet and Sulem [7]. We have run for longer times than 
the results shown in [7], Plots for this run are given in figures 9 - 10 . 

Run IV: Random forcing. 

In this run we consider the problem when the equations are forced in a range of 
low Fourier modes. For the forced problem there appears to be no easy way to obtain 
a sharp bound on the maximum of the vorticity. We have found experimentally that 
when the forcing is chosen to be 0(1) the solution grows and does not remain 0(1). 
For example in figure 11 we show the results of a run in which the forcing f is 
chosen so that l/loo = 1 and in which the initial vorticity is Eero. In particular the 
amplitudes and phases of the the fourier components of the forcing were chosen as 


'Cl Cl 

Ci«i 

ClC 2 

ci< 2 

ClC S 

Cl <3 " 


r+12 

-4 

-16 

+4 

+20 

- 12 - 

«lCl 

< 1*1 

<ic 2 

< 1*2 

<lCs 

< 1*3 


+8 

-24 

-20 

+12 

+28 

+4 

CjCx 

C 2 *l 

C 2 C 2 

C 2 S 2 

C 2 C 3 

C 2*3 

— c 

+12 

-8 

+32 

+24 

+8 

+36 

« 2 Ci 

<a*i 

<JC 2 

< 2*2 

»2C3 

<2*3 


-12 

+12 

+4 

+32 

-4 

-16 

C 3 C 1 

Cs*l 

C 3 CJ 

Cs< 2 

C 3 C 3 

Cs*s 


-4 

-36 

-16 

-36 

-24 

+4 

■*SCi 

»S*1 

<SC 2 

<3*2 

<SC3 

< 3«3 - 


.-20 

-20 

-36 

+8 

-28 

+12. 
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where the scaling factor C t was chosen to ensure that \f\ 00 = 1. In this run 
v = 10" 8 , Ni= N 2 = 256 and (CFL^, CFL opl , CFL mM ) = (.5, .8,1.). In figure 
11 b we have made plots of the maximum norms 

MAX(U) = max(|it(-,.,t)| 00 ,|t>(.,-,<)(oo) 

and 

MAX(W) = |«( v ,0|oo. 

Even by time t = 200 the solution continues to grow. 

In contrast when the forcing is chosen to be 0( v) we do not see growth in the 
|u>(-, *,£)|oo- This observation is presented in figures 12 - 13 where we have made 
runs with v = .5 x 10 -8 (JV = 128) and v = .5 x 10 -4 ( N = 256). The initial 
conditions for |u>| were defined by the matrix of coefficients given above but in this 
case the constant C, was chosen so that !<*>(•, •, 0)|oo = 1. The forcing was constant 
in time and defined from the relation 

i/Aw + f = 0. 


4 Discussion 

We have shown that for both two and three dimensional flows the minimum 
scale, A m j„, is essentially proportional to */ 1 / a /|Du| 00 . We now relate this 

minimum scale to the decay rate of the energy spectrum. 

Let us assume that at a given time t the energy spectrum has a power law 
behavior in some range of wave numbers, the inertial range , 

E(k) ~ k~™, kt = 0(1) < k < 0(1/A m , n ). 


We have proved that the quantities 
i '"- 1 i /”- 1 


\Du\ ( 


= 


UH 


p + 1 
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remain of order one, provided they are initially so. When this order one bound is 
achieved we call the flow maximal digtipaiive. Note that 2 vH\ is the rate of energy 
dissipation, e, and that the rate of enstrophy dissipation rj is bounded by 

Assuming that the leading order contribution to the integral for H*(t) is deter- 
mined from wave numbers in the inertial range and using the power law behavior 
for E{k) we obtain 


yP~ l 


|Du| 


P+1 P 


Hl{t) 




k 2, E ( k ) dk „ 


\Du\ 


r+i 


- 0 - 1/2 


0 + 8/2 


In two space dimensions we know that IDu^ is essentially bounded by the 
maximum norm of the initial vorticity. Let us thus assume that the initial values 


are scaled so that IDu]^ is of order one. In this case 

I /*" 1 


MH 


P + 1 


fli(0 - 


(4.1) 


and for maximal dissipative flows it follows that E(k) ~ k~ a , the power law behavior 
predicted by the Batchelor-Kraichnan theory [3] [16]. 

In three dimensions if we speculate that [Du^ = 0(i/~ y ), then 

l/P ~ 1 ff ! ^ tt(0+1/2)+0-3/2 

P«L' +1 ' 

In this way for maximal dissipative flows, we obtain a relation between the power 
law behaviour of the energy spectrum and the size of 


0 = 


3-7 
2 + 27' 


When (Du)^ = 0(i /~ l ^ 2 ), and 7 = 1/2, we obtain 0 = 5/6 and E(k) = A -5 / 8 , the 
power law behavior predicted by Kolmogoroff [15]. 

We now return to the numerical results of the random initial data runs, refer- 
ence figures 2-4 and 7 - 8 . The initial conditions were chosen so that E(k) ~ k~ s . 
The numerical results show that this A -8 power law seems to remain over an initial 
time interval. The numerical results further indicate that as the flow evolves, the 


- 32 - 


quantities v p 3 which behave like v p H p +\, slowly decrase and the energy spectrum 
steepens. In a later regime the flow is dominated by the presence of large regions 
of relatively constant vorticity. There seems to be some evidence from the figures 
to suggest that the quantities 3* decay by a factor of about A m ,- n = t/ 1 / 2 , and that 
i(k) ~ A" 15 , E(k) ~ Jfe" 4 . When 


.p-i 


|Cu| 


P+1 


H 2 p (t) 


1 1 


1/2 


the argument in (4.1) predicts /3 = 2, and is thus consistent with the numerical 
results. These results are in agreement with Saffman’s theory for two dimensional 
turbulence [22], 

In figure 1 we outline an hypothesized behaviour of i/ p J 2 (<) for the decay of 
two-dimensional turbulence. In the first stage of development of the flow, i / p J* may 
show an overall increase as the flow evolves to a state of maximal dissipation. (Of 
course, depending on the initial data, this maximal dissipative state may never be 
reached.) This dissipation rate can not continue for a long time interval but must 
decrease. The power law then slowly changes from k~ s to a more rapid decay. The 
flow becomes organized into coherent structures, a regime with i / p J p " 1/2 (?) 
and where Saffman’s theory would predict E{k) ~ At -4 . This regime presumably 
exists for long times, since the viscosity now plays a minimal role. This scenario is 
suggested by our computations and other similar ones. In particular, Brachet and 
Sulem [7] show high resolution computation with initial data 


E{k) 


similar to the one presented in this paper (reference figures 9 - 10 ). They found an 
increase in the energy power law reaching a maximum at about At -9 . At this stage 
the rate of enstrophy dissipation is maximum, in accordance with our analytical 
results. 

It is conceivable that a similar scenario is present in the decay of three dimen- 
sional turbulence. In two dimensions when large coherent structures are formed 
the main contribution to H*(t) comes from one dimensional layers separating these 
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structures. When the solutions across the layers have a simple structure, one can 
argue that that 

i /*- 1 , i / 1 / 2 

FFT^fW ~ ■'min ~ ==— 1 7, ■ 

l^»L 

Correspondingly, in three dimensions the same argument can be made assuming that 
the regions of rapid variation are concentrated along two dimensional structures of 
width Xrnin- Thus in either case we obtain a new relation between the power law 
behaviour of the energy spectrum and the size of jDul^ = 0(t/ _ ' r ): 


\Du\ 


In two dimensions 7 = 0 and we again obtain E{k) ~ k~*. If, in three dimensions 
l-Dul^ ~ i/ 1 / 2 , we obtain E(k ) ~ fc -8 / 3 . Large three dimensional simulations are 
necessary to confirm the validity of the assumption made on the size of IDu^, on 
the time evolution of H*[t) and the sharpness of our estimates. 



-34- 


References 

[1] J. Barker, A Numerical Experiment on the Structure of Two-Dimensional Tur- 
bulent Flow, Ph.D. Thesis, Applied Math. Dept., Caltech (1982). 

[2] C. Basdevant, Technical Improvements for Direct Numerical Simulation of Ho- 
mogeneous Three-Dimensional Turbulence , J. Comp. Physics 50 (1983), pp. 
209-214. 

[3] G.K. Batchelor, Computation of the Energy Spectrum in Homogeneous Two- 
Dimensional Turbulence, Phys. Fluid. Snppl. II (1969), pp. 233-239. 

[4] R. Benzi, G. P&ladini, S. Patamello, P. S ant an gel o, A. Valpiani, Intermittency 
and Coherent Structures in Two-Dimensional Turbulence, J. Phys. A: Math. 
Gen. 19 (1986), pp. 3771-3784. 

[5] M.E. Brachet, D.I. Meiron, S.A. Orszag, B.G. Nickel, R.H. Morf, U. Frisch, 
Small-scale Structure of the Taylor-Green Vortex, J.F.M. 130 (1983), pp. 411- 
452. 

[6] M.E. Brachet, M. Meneguzzi, P.L. Salem, Small Scale Dynamics of High Reynolds 
Number Two Dimensional Turbulence, snbm. to Phys. Rev. Let. 

[7] M.E. Brachet, P.L. Salem, Direct Numerical Simulation of Two Dimensional 
Turbulence, Proc. 4th Beer Shevaof MHD Flows and Turbulence, March 1984, 
Israel. 

[8] B. Fornberg, A Numerical Study of i-D Turbulence, J. Comp. Phys., 25 (1977), 
pp. 1-31. 

[9] D.G. Fox, S.A. Oiszag, Pseudospectrol Approximation to Two Dimensional 
Turbulence, J. Comp. Phys., 11 (1973), pp. 612-619. 

[10] D. Gilbarg, N.S. Tkndinger, Elliptic Partial Differential Equations of Second 
Order, Springer Verlag, (1983). 



- 35 - 


[11] J.R. Herring, J.C. McWilliams, Comparison of Direct Simulation of Two Di- 
mensional Turbulence with Two-Point Closure: the Effects of Intermittency, J. 
Fluid Mech., 153 (1985), pp. 229-242. 

[12] J.R. Herring, S.A. Orssag, R.H.Kraichnan, D.G. Fox, Decay of Two-Dimensional 
Turbulence, J. Fluid Mech., 66 (1974), pp, 417-444. 

[13] L. Hormander, The Boundary Problems of Physical Geodesy, Arch. Rational 
Mech. Anal., 62 (1976), pp. 1-52. 

[14] M.Y. Hussaini, T.A. Zang, Spectral Methods in Fluid Dynamics, Ann. Rev. 
Fluid Mech., 10 (1987), pp. 339-367. 

[15] A.N. KolmogorofF, C.R.Acad.Sci. U.S.S.R. 30 (1941), pp. 301. 

[16] R. Kraichnan, Inertial Ranges in Two Dimensional Turbulence, Phys. Fluids, 
10 (1967), pp. 1417-1423. 

[17] H.O. Kreiss, W.D. Henshaw, L.G. Reyna, On the Smallest Scale for the Incom- 
pressible Navier-Stokes Equations in Two Space Dimensions, IBM Research 
Report No. 13204 (1987). 

[18] H.O. Kreiss, J. Oliger, Comparison of Accurate Methods for the Integration of 
Hyperbolic Equations, Tellus 24 (1972), pp. 199-215. 

[19] D.K. Lilly, Numerical Simulation of Developing and Decaying of Two Dimen- 
sional Turbulence, J. Fluid Mech., 45 (1971), pp. 395-415. 

[20] S.A. Orssag, Analytical Theories of Turbulence, J. Fluid Mech., 41 (1970), pp. 
363-386. 

[21] S.A. Orseag, Numerical Simulation of Incompressible Flows Within Simple 
Boundaries: I. Galerkin (Spectral) Representation, Stud. Appl. Math., 50 
(1971), pp. 293-327. 

[22] P.G. Saffman, On the Spectrum and Decay of Random Two Dimensional Vortic - 
ity Distributions at Large Reynolds Numbers, Studies Appl. Math., 50 (1971), 
pp. 377-383. 


- 36 - 


Vorticity t= 0.00 v=0.1000E-03 Vorticity i= 20 .00 v=0 .1000E-03 
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Figure 2 Contour plots of the vorticity for random initial data, u = 10 4 , N = 256. 
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Vorticity Maz = 0.114 t= 0 .00 v=0 .1000E-03 Vorticity Maz = 0. 114 t= 20.02 u=0 .1000E-03 



Vorticity Max = 0.116 t= 50.01 v=0. 1000E-02 Vorticity Maz = 0 .116 MOO. 01 v=0.1000E-03 



Figure 3 Surface plots of vorticity spectrum for random initial data, u = lO -4 , A" = 256. 


Figure 4 Results for random initial data, v = 10 


N = 256. 
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Vorlicily t= 0.00 v=0.1000E-03 Vorlicity t= 20 .02 v=0 .1000E-03 
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Figure 5 Contour plots of the vorlicity for random initial data, u - 10 4 , N = 128. 
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Figure 6 Results for random initial data, = 10 4 , A r = 128. 


ORIGINAL PAGE IS 
OF POOR QUA 1.1733 





-41- 

ORIGINAL PAGE IS 
QE POOR QUALITY 


Korf icily 1= 0 .00 v-0 .1000E-04 Voriicitv t= SO .00 v=0 1000E-04 



Figure 7 Contour plots of the vorticity for random initial data, v = 10“*, N = 512. 
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Figure 8 Results for random initial data, v = 10 5 , N = 512. 
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Figure 9 Contour plots of the vorticity for smooth random initial data, i/ = 2xlO -s ,JV = 256. 
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Figure 10 Results for smooth random initial data, 1 / — 2 x 10 5 , N = 256. 
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Figure 12 Plots of norms for random forcing, \f\oc = 0(i/), 

v = .5 x 1CT 3 , N = 128 and i/=.5x 10" 4 , N = 256. 
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Figure 13 Contour plots of vorticity spectrum for random forcing, !/!* = O(u), 

„ = .5 x 10“®, N = 128 and i/ = .5x 10' 4 , N = 256 
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Figure 14 Hypothetical behaviour for v r J*(t) 
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